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PREFACE 


This  report  outlines  recent  efforts  to  use  finite  element  techniques  for  solving  the  wave  equation 
of  underwater  acoustic  propagation  in  large  ocean  regimes.  This  effort  is  a  part  of  a  computationally 
intensive  probabilistic  acoustics  program  whose  major  current  goal  is  to  develop  models  for  the  propaga¬ 
tion  of  the  moments  of  the  acoustic  field  in  regions  where  the  boundaries  are  random  surfaces.  With 
the  rapidly  increasing  computer  power  of  large  main-frame  computers,  indeed  even  large  minicomput¬ 
ers  augmented  with  powerful  array  processors,  the  possibility  of  using  an  exact  technique  as  the  finite 
element  method  (FEM)  for  solving  the  wave  equation,  even  in  such  complicated  regimes,  needs  to  be 
addressed.  In  the  near  term  this  method  is  expected  to  be  especially  applicable  for  benchmark  or 
exact-solution  calculations  subsequently  used  to  evaluate  large  complex  computer  codes,  since  the 
length  of  time  that  such  calculations  take  is  not,  in  general,  the  governing  parameter.  Future  advances 
in  computer  design  can  make  the  FEM  also  attractive  for  solving  problems  in  realistic  scenarios  and  for 
use  in  a  production  mode.  For  relatively  small  source-to-receiver  ranges  at  low  frequency,  the  FEM 
can  already  be  used  as  a  predictor  of  transmission  loss  in  a  complicated  environment. 
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THE  NUMERICAL  SOLUTION  OF 
UNDERWATER  ACOUSTIC  PROPAGATION  PROBLEMS  USING 
FINITE  DIFFERENCE  AND  FINITE  ELEMENT  METHODS 


INTRODUCTION 

In  this  report,  we  discuss  various  aspects  of  the  numerical  solution  of  underwater  acoustic-wave 
propagation  problems.  The  propagation  of  acoustic  energy  in  the  ocean  involves  the  interaction 
between  acoustic-wave  propagation  in  fluids  and  stress-wave  propagation  in  solids.  Thus,  a  general 
mathematical  model  involves  the  coupling  of  the  acoustic-wave  equation  with  the  elastic-wave  equation 
and  the  specification  of  suitable  interface  and  boundary  conditions.  Only  simple  wave-propagation  prob¬ 
lems  can  be  solved  analytically.  Hence,  the  approximate  solution  of  time-harmonic  and  time-dependent 
models  in  two  and  three  dimensions  is  important  to  treat  effectively  acoustic  propagation  in  a  general 
ocean  environment.  Note  that  we  are  considering  linear,  forward,  deterministic  propagation  problems 
here.  However,  many  of  the  methods  may  also  be  applicable  to  nonlinear,  inverse,  and  stochastic  prob¬ 
lems. 


Various  computational  approaches  have  been  developed  and  applied  to  simplified  propagation 
models.  These  include  parabolic-equation  and  normal-mode  models,  asymptotic  methods,  and  others; 
see,  for  example,  Ref.  1  for  a  survey  of  various  models  and  numerical  techniques.  Although  each  of 
these  techniques  can  be  quite  effective  under  suitable  assumptions,  there  are  many  important  problems 
for  which  it  is  necessary  to  treat  the  complete  wave-propagation  model  described  above  in  the  low  to 
intermediate  frequency  range.  Such  models  can  include  lateral  inhomogeneities,  multiple  irregular 
interfaces  and  boundaries,  full  angle  propagation,  and  backscattering.  This  occurs,  for  example,  when 
the  ocean  bottom  must  be  taken  into  account,  such  as  in  shallow-water  propagation  and  in  deep  water 
at  very  low  frequencies.  The  interaction  of  acoustic  and  seismic  waves  with  a  complicated  ocean  bot¬ 
tom  is  an  important  and  difficult  problem.  Another  important  example  of  a  complicated  propagation 
problem  occurs  when  a  layer  of  ice  is  present  on  the  ocean  surface. 

Finite  difference  and  finite  element  methods  have  proved  to  be  very  effective  techniques  for  solv¬ 
ing  approximately  boundary  and  initial-value  problems  of  the  type  described  above.  However,  to  prop¬ 
erly  resolve  the  waves  it  is  necessary  to  decrease  the  spatial  mesh  sizes  as  the  frequency  increases.  This 
can  result  in  problems  with  very  large  numbers  of  degrees  of  freedom  as  the  frequency  and/or  spatial 
dimensions  increase.  The  size  of  the  problems  that  can  be  effectively  treated  numerically  depends  on 
various  factors,  such  as  the  computer  power,  as  well  as  the  mathematical  and  modeling  techniques 
available.  We  shall  discuss  these  and  other  issues  in  more  detail  later  and  consider  various  possibilities 
for  treating  large,  complicated  ocean-propagation  problems  more  efficiently. 

We  close  this  section  by  outlining  the  remainder  of  this  report.  In  the  second  section,  we  con¬ 
sider  a  model  problem  involving  the  two-dimensional  Helmholtz  equation  with  a  variable  sound  speed; 
this  was  treated  using  a  finite-element  code  implemented  at  NRL  on  the  VAX  11/780  and  modified  to 
treat  underwater  acoustic-propagation  problems.  We  briefly  describe  some  features  of  the  finite- 
element  algorithm  and  the  treatment  of  radiation  boundary  conditions.  A  distinctive  feature  of  the 
code  is  the  implementation  of  a  recently  developed  iterative  method  [2]  for  solving  the  resulting  large, 
sparse,  indefinite,  non-self-adjoint  system  of  equations.  This  allowed  for  the  efficient  solution  of  over 
35,000  complex  equations  on  a  relatively  small  computer,  since  large  matrices  did  not  have  to  be  stored 
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or  inverted.  We  describe  some  of  the  numerical  results  obtained  after  applying  this  code  to  the  model 
problem.  Furthermore,  we  discuss  additional  modifications  that  can  be  made  to  the  code  to  improve  its 
efficiency  and  extend  its  applicability  to  more  general  propagation  models. 

In  the  third  section,  we  consider  the  general  situation  of  the  coupled  acoustic/elastic  wave  equa¬ 
tion  in  two  and  three  dimensions.  We  discuss  finite-difference  and  finite-element  methods  for  solving 
both  the  time-harmonic  and  time-dependent  model.  The  numerical  techniques  for  solving  the  time- 
harmonic  and  time-dependent  problems  are  very  different,  as  are  the  numerical  difficulties  that  are 
encountered.  Various  issues  are  considered  that  are  important  in  determining  the  size  of  the  problem 
that  can  be  adequately  treated.  Finally,  in  the  fourth  section  we  summarize  our  findings. 

A  MODEL  PROBLEM 


In  this  section  we  describe  a  model  underwater  acoustic  propagation  problem  based  on  the  two- 
dimensional  Helmholtz  equation  with  a  variable  sound  speed.  We  also  describe  results  obtained  after 
implementing  a  finite-element  code  and  modifying  it  to  treat  this  propagation  model.  Finally,  we  indi¬ 
cate  several  ways  of  improving  the  capabilities  of  this  code. 


The  Model  and  the  Numerical  Algorithm 


The  Helmholtz  equation  for  a  cylindrically  symmetric  geometry  and  a  harmonic  source  is  given  by 

A  u(r,z)  +  Kq  n2(r,z)u(r,z )  =  0, 


where  z  is  the  depth  measured  downard  from  the  ocean  surface,  r  is  the  range,  u  ( r,z )  is  the  acoustic 
pressure, 


a  1  s 

Am  =  —  — 
r  dr 


9m 

'dr 


is  the  Laplacian  in  cylindrical  coordinates,  the  reference  wave  number  is  K0  =  lirf/co,  f  is  the  source 
frequency,  c0  is  a  reference  sound  speed,  the  refractive  index  n(r,z)  =  cjc{r,z),  and  c(r,z)  is  the 
sound  speed.  For  simplicity,  we  consider  a  region  of  the  ocean  with  a  flat  surface  and  a  flat  bottom,  so 
that  the  region  D  is  a  rectangle  (see  Fig.  1).  We  assume  an  ocean  depth  of  5000  m,  an  initial  range 
given  by  r  —  R 1,  and  a  final  range  given  by  r  =  R2. 


Our  boundary  conditions  are  given  by 

Top:  u(r,0)  =  0,  (la) 

Bottom:  ^-(r,  5000)  =  0,  (lb) 
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and 


Left  side:  u{R\,z )  =  g(z),  (lc) 

where  g(z)  is  a  specified  initial  pressure  field.  The  boundary  condition  at  r  =  R2  is  chosen  so  as  to 
model  the  outgoing  radiation  of  energy.  Since  there  is  no  backscattering,  we  may  model  this  outgoing 
radiation  condition  by  choosing  a  suitable  artificial  dissipation,  id(r),  in  the  region  RDIS  <  r  <  R2, 
where  dir)  is  a  real-valued  nonnegative  function  to  be  specified  later.  The  right-hand  boundary  condi¬ 
tion  is  now  defined  by 

Right  side:  uiR2,z)  =  0.  (Id) 

We  shall  see  later  that  there  are  various  alternative  methods  for  modeling  the  outgoing  radiation  condi¬ 
tion. 


Our  goal  is  to  solve  approximately  the  following  elliptic  partial  differential  equation  in  the  region 
D,  combined  with  boundary  conditions  (la)  through  (Id): 

A  u(r,z)  +  Kon2(r,z)[l  +  id(r)]u(r,z )  =  0.  (2) 

We  have  run  our  computer  code  with  the  following  choices  of  parameters  and  functions: 

V^o 


g(z)  = 


exp  |  -  K$  (z-z0)2/4j  -  exp  -  K$  (z+z0)2/4| 


Cq  =  1500  m/s,  source  depth  z0  —  2500  m, 


*  =  —  a2(z—zi)2,  a  =  10  7  m2/s. 


c2(z) 


depth  of  sound  speed  minimum  z\  —  2500  m,  C\  =  1500  m/s,  i?l  =  1  m,  and  R 2  has  been  chosen 
thus  far  between  3000  and  25,000  m.  The  frequency  / has  been  chosen  thus  far  between  3  and  10  Hz. 
This  particular  index  of  refraction  causes  a  "focus”  at  about  20  km.  Note  that  while  the  initial  Gaussian 
pressure  field  g(z )  does  not  exactly  satisfy  boundary  conditions  (la)  and  (lb),  these  boundary  condi¬ 
tions  are  sufficiently  closely  satisfied  for  our  choice  of  parameters  that  this  causes  no  computational 
difficulties.  This  model  was  previously  run  using  a  parabolic  equation  method  implemented  at  NRL 
[3],  although  the  bottom-boundary  condition  in  Ref.  3  differs  from  Eq.  (lb).  We  have  found  empiri¬ 
cally  that  a  convenient  functional  form  for  the  dissipation  term,  idir),  is  given  by 


dir) 


ef3(r  —  RDIS) — !  for  RDIS  <r<R2, 
0,  for  r  <  RDIS, 


with  (3  and  RDIS  suitably  chosen,  although  other  simple  functions  can  work  about  as  well.  The  pur¬ 
pose  of  this  dissipation  term  is  to  attenuate  the  wave  while  at  the  same  time  minimizing  the  reflection 
due  to  the  dissipation.  Note:  There  is  no  difficulty  in  treating  range-  and  depth-dependent  sound 
speeds.  Also,  R  1  can  be  chosen  quite  large.  The  input  data,  g ,  may  be  specified  as  the  result  of  a 
(long-range)  parabolic  equation  run  or  some  other  numerical  or  asymptotic  method.  Finally,  the  model 
problem  may  be  readily  generalized  to  include  complicated  geometries,  boundary  conditions,  and  inter¬ 
faces.  This  latter  point  is  discussed  in  more  detail  later. 


The  numerical  algorithm  we  have  employed  to  solve  approximately  the  boundary-value  problem 
given  by  Eqs.  (1)  and  (2)  is  based  on  the  finite  element  method.  We  shall  not  go  into  a  technical  discus¬ 
sion  of  the  finite-element  method,  since  it  is  described  extensively  in  the  literature  (see,  e.g.,  Ref.  4). 
We  merely  point  out  that  the  finite  element  method  is  based  on  replacing  the  given  boundary-value 
problem  by  an  equivalent  variational  problem  and  then  approximating  the  variational  problem  by  use  of 
a  convenient  finite-dimensional  space  of  functions.  Typically,  this  space  of  functions  consists  of 
sufficiently  smooth  piecewise  polynomials  defined  with  respect  to  a  partitioning  of  the  computational 
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domain  into  simple  subsets  called  elements.  This  reduces  the  variational  problem  to  that  of  solving  a 
finite  number  of  linear  equations.  As  the  diameter  of  the  largest  element  decreases,  the  approximate 
solution  converges  to  the  exact  solution,  but  the  number  of  equations  to  be  solved  increases.  The  com¬ 
puter  code  we  have  implemented  to  solve  the  problem  is  based  on  continuous  piecewise  linear  func¬ 
tions  defined  on  right  triangles.  The  code  can  use  either  uniform  or  variable  mesh  sizes.  Let  us 
observe  that  finite  difference  methods  can  be  used  just  as  effectively  to  solve  the  problem  approxi¬ 
mately.  Finite  element  methods,  however,  are  more  effective  for  treating  complicated  boundaries, 
interfaces,  and  boundary  conditions. 

An  important,  distinctive  feature  of  our  numerical  algorithm  and  computer  code  is  the  incorpora¬ 
tion  of  a  recently  developed  iterative  method  [2]  for  solving  the  resulting  sparse  system  of  linear  equa¬ 
tions.  The  solution  of  this  system  of  equations  is  the  most  expensive  part  of  the  computation.  It  is 
well-known  that  iterative  methods  are  in  general  considerably  more  efficient  than  direct  methods  (i.e., 
those  based  on  Gaussian  elimination)  for  large  problems,  with  respect  to  both  storage  and  computa¬ 
tional  speed.  However,  iterative  methods  have  typically  been  developed  and  analyzed  for  positive 
definite,  symmetric  problems  (see,  e.g.,  Ref.  5).  Neither  of  these  properties  holds  for  the  problems 
currently  under  investigation.  An  iterative  method  based  on  the  preconditioned  conjugate  gradient 
method  was  described  in  Ref.  2  for  a  class  of  problems  including  the  time-harmonic  problems  discussed 
in  this  report.  The  preconditioner  is  based  on  one  sweep  of  symmetric  successive  overrelaxation 
(SSOR),  although  other  preconditioners  are  being  investigated.  Its  implementation  has  resulted  in  a 
dramatic  increase  in  storage  capabilities  and  a  dramatic  decrease  in  computer  time  with  respect  to  direct 
solvers,  since  large,  sparse  matrices  do  not  have  to  be  stored  or  inverted* 

Numerical  Results 

We  next  briefly  describe  results  obtained  after  applying  the  finite  element  code  to  the  boundary- 
value  problem  given  by  Eqs.  (1)  and  (2).  This  work  concentrated  on  the  following: 

•  Implementing  a  finite  element  code  developed  elsewhere  on  the  VAX  11/780  in  the  Large 
Aperture  Acoustics  Branch  at  the  Naval  Research  Laboratory  and  modifying  the  code  to 
treat  propagation  problems  such  as  the  aforementioned  one; 

•  Improving  the  code  so  as  to  increase  its  speed  and,  particularly,  its  storage  capabilities; 
and 

•  Testing  important  quantities,  such  as  the  range  of  artificial  dissipation  and  the  number  of 
grid  points/wavelength  needed  for  a  prescribed  accuracy,  as  well  as  the  CPU  time  for  vari¬ 
ous  choices  of  parameters. 

It  is  important  to  emphasize  that  we  are  solving  large  problems  on  a  relatively  small  computer. 
We  are  now  able  to  store  and  solve  over  35,000  complex  equations  on  the  VAX  11/780.  (When  the 
code  was  originally  implemented,  we  were  limited  to  about  12,000  equations.)  Furthermore,  the 
subroutines  used  for  the  iterative  method  have  been  improved  so  as  to  make  them  more  efficient  for 
use  on  a  vector  computer  or  an  array  processor.  There  is  still  a  great  deal  that  can  be  done  to  increase 
the  efficiency  of  the  code.  This  is  outlined  later. 

To  determine  the  length  of  artificial  dissipation  and  the  number  of  grid  points/ wavelength  needed, 
we  used  measures  based  on  the  average  volume  intensity  and  average  line  intensity  of  the  computed 
solution.  We  can  plot  the  solution  as  a  surface,  and  we  can  plot  the  transmission  loss  vs  range.  We 
have  determined  empirically,  using  these  measures  and  plots,  that  when  the  length  of  the  dissipation 
layer  is  about  two  to  four  wavelengths,  there  is  no  significant  deterioration  in  the  solution.  As  for  the 
resolution  of  the  waves,  we  observed  that  a  minimum  of  8  grid  points/wavelength  is  necessary  to 
obtain  a  meaningful  solution  when  a  uniform  grid  size  is  used. 
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To  get  an  idea  of  the  CPU  time  for  typical  runs,  we  consider  the  following  sample  test  runs  at  3, 
5,  and  10  Hz  for  R 2  =  2001  m,  using  approximately  8  and  12  grid  points/wavelength  (see  Table  1). 
Note  that  as  the  frequency,  /,  is  multiplied  by  a  factor  C,  the  number  of  equations  is  multiplied  by 
approximately  C2  (to  resolve  the  waves  in  two  dimensions).  Furthermore,  the  number  of  iterations 
generally  increases  only  slightly.  This  gives  a  rough  idea  of  how  the  CPU  time  increases  with  frequency 
for  fixed  range.  We  also  observe  that  as  the  mesh  size  decreases  from  a  coarse  mesh  (8 
points/wavelength)  to  a  finer  mesh  (12  points/ wavelength),  the  number  of  iterations  increases  very 
slightly  or  can  even  decrease  even  though  the  number  of  equations  increases.  This  occurs  because  the 
system  of  equations  is  better  conditioned  as  the  discrete  model  better  approximates  the  continuous 
model.  Further  study  is  needed  to  assess  the  accuracy  of  the  discrete  model  with  respect  to  changes  in 
all  the  parameters. 


Table  1 


Case 

/ 

(Hz) 

Points/ 

Wavelength 

Number  of 
Equations 

Number  of 
Iterations 

CPU  Time 
(min) 

(a) 

3 

8 

2,673 

169 

8 

(b) 

5 

8 

7,809 

287 

39 

(c) 

10 

8 

28,944 

287 

150 

(d) 

3 

12 

5,929 

196 

21 

(e) 

5 

12 

17,425 

216 

84 

Next,  suppose  that  the  frequency  and  grid  points/wavelength  are  fixed  but  the  range  increases  by 
a  factor  C.  Hence,  the  time  for  each  iteration  increases  by  a  factor  C.  Furthermore,  the  number  of 
iterations  will  increase  typically  by  a  factor  C',  where  1  <  C'  <  C.  For  example,  consider  case  (a)  in 
Table  1  (3  Hz  over  2000  m).  We  ran  it  over  25,000  m,  so  that  N  -  32,481.  The  CPU  time  was  414 
min.  (The  number  of  iterations  was  multiplied  by  a  factor  of  4.)  In  terms  of  storage  on  the  VAX,  the 
current  version  of  the  code  can  treat  5  Hz  over  10,000  m  or  10  Hz  over  2500  m  using  a  uniform  mesh. 

We  next  describe  a  factorization-mesh  grading  procedure  developed  to  treat  longer  ranges  without 
increasing  the  storage  capabilities.  This  method  consists  first  of  factoring  out  e  °r  from  the  solution. 
This  results  in  a  smoother  solution  as  r  increases  and  thus  allows  for  longer  range  step  sizes  as  we 
proceed  away  from  the  origin.  This  is  consistent  with  the  approach  taken  when  the  parabolic  equation 
method  is  employed,  where  it  is  typically  observed  that  range  steps  several  wavelengths  long  suffice  for 
accurate  far-field  solutions.  This  approach  is  also  analogous  to  a  method  developed  and  analyzed  in 
Ref.  6  in  connection  with  the  Helmholtz  equation  exterior  to  a  bounded  obstacle.  In  the  current  code, 
we  implemented  this  factorization-mesh  grading  procedure  using  larger  range  steps  with  increasing  r. 
Due  to  lack  of  time,  we  were  unable  to  study  its  effectiveness  comprehensively.  However,  preliminary 
results  indicated  that  we  could  substantially  reduce  the  number  of  grid  points  without  losing  accuracy. 
For  example,  we  ran  a  problem  with  3  Hz  over  30,000  m  using  this  procedure  and  were  able  to  reduce 
the  number  of  range  points  by  nearly  a  factor  of  three  compared  to  a  uniform  grid  with  8 
points/ wavelength. 

Suggested  Improvements  and  Extensions 

We  conclude  this  section  by  outlining  several  modifications  that  can  be  made  to  increase  the  capa¬ 
bilities  of  the  present  code  and  shed  more  light  on  its  efficiency  and  generality. 

(1)  There  is  no  difficulty  in  changing  the  initial  pressure  field,  g(z),  and  sound-speed  profile, 
c(r,z),  in  the  computer  code.  Hence  the  code  could  be  readily  employed  to  solve  a  propagation  prob¬ 
lem  for  which  the  exact  solution  is  known.  For  example,  g(z)  could  be  generated  by  use  of  a  normal¬ 
mode  solution.  This  could  be  useful  in  obtaining  more  detailed  information  regarding  the  accuracy  of 
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The  two  vertical  boundary  conditions  can  also  be  modeled  in  different  ways.  For  the  left  bound¬ 
ary  condition  we  might  assume,  as  in  the  previous  section,  that  the  field  is  specified,  perhaps  as  the 
result  of  a  long-range  run  with  a  simplified  propagation  model  used.  Alternatively,  the  loading  might 
be  such  that  a  plane  of  symmetry  exists.  In  this  case  the  particle  motions  are  zero  normal  to  the  plane 
of  symmetry  on  this  boundary.  The  right  vertical  boundary  requires  the  imposition  of  an  outgoing- 
radiation  (nonreflecting)  boundary  condition.  This  is  an  active  area  of  research  and  can  have  an  impor¬ 
tant  impact  on  the  computations.  We  shall  consider  alternative  formulations  of  this  radiation  condition 
later. 


The  resulting  computational  model  can  be  made  discrete  by  the  use  of  either  finite  difference  or 
finite  element  techniques.  This  typically  results  in  a  large,  sparse  system  of  linear  equations.  The  prob¬ 
lem  size  depends  mainly  on  the  frequency  and  the  size  of  the  computational  domain.  If  we  assume  for 
simplicity  a  uniform  grid  size,  it  is  clear  that  if  we  multiply  the  length  of  the  domain  by  m,  then  the 
number  of  equations  is  multiplied  by  m.  However,  multiplying  the  frequency  by  m  results  in  multiply¬ 
ing  the  number  of  equations  by  at  least  m 2  and  m3  in  two  and  three  dimensions,  respectively,  since  the 
waves  must  be  resolved  in  each  spatial  direction.  The  frequencies  and  dimensions  of  interest  depend, 
of  course,  on  the  particular  physical  problem. 

We  thus  see  that  the  problem  size  can  be  quite  large,  particularly  as  the  frequency  increases. 
(Asymptotic  methods,  such  as  ray  tracing,  can  be  effective  for  high  frequencies.  However,  it  is  not 
clear  in  general  for  which  frequencies  these  methods  yield  reliable  results.)  In  the  remainder  of  this 
section,  we  discuss  various  issues  relating  to  the  size  of  the  propagation  models  that  might  reasonably 
be  solved  numerically.  We  first  consider  the  computer  power  that  is  currently  available  and  that  is 
expected  to  become  available  in  the  next  few  years.  We  then  discuss  some  mathematical  techniques 
tha*  can  have  an  important  bearing  on  the  solution  of  the  computational  models.  The  mathematical 
methods  for  treating  time-harmonic  and  transient  models  are  very  different  and  will  be  discussed 
separately. 

Computational  Power 

Because  of  the  increasing  need  for  large-scale,  high-speed  computers  in  many  different  areas  of 
science  and  engineering,  there  is  a  great  deal  of  work  proceeding  with  the  aim  of  greatly  increasing 
computing  power.  We  first  briefly  discuss  the  increased  computer  power  obtainable  using  vector 
machines  such  as  the  CRAY-1  or  CYBER  205.  We  observe  in  this  connection  that  the  finite-element 
algorithm  we  described  previously  is,  for  the  most  part,  vectorizable.  The  main  exception  is  the  SSOR 
preconditioner  used  in  the  iterative  method.  We  discuss  later  alternative  preconditioners  whose  imple¬ 
mentation  can  make  the  iterative  method  completely  vectorizable.  We  will  attempt  to  compare  the 
computational  power  of  the  CRAY-1,  CRAY  X-MP,  and  CRAY-2  with  that  of  the  VAX  11/780  used 
for  the  computations  discussed  in  the  previous  section.  These  comparisons  are  based  on  the  best 
knowledge  and  conjectures  available  to  the  author  at  this  time. 

The  CRAY-1,  currently  in  use,  has  a  memory  of  between  2  x  106  and  4  x  106  64-bit  words.  The 
memory  available  to  us  on  the  VAX  was  about  1  x  106  32-bit  words.  Note  that  the  use  of  64  bits  on 
the  CRAY  is  roughly  equivalent  to  the  use  of  double  precision  on  the  VAX.  The  CRAY-1  is  approxi¬ 
mately  40  times  as  fast  as  the  VAX.  The  CRAY  X-MP  is  scheduled  to  be  operating  sometime  in  1984 
and  is  to  be  about  twice  as  fast  as  the  CRAY-1.  Finally,  the  CRAY-2  is  to  be  a  substantial 
improvement  over  the  CRAY-1  in  both  memory  and  speed.  Specific  details  at  this  point  are  not 
known.  However,  available  information  indicates  an  improvement  in  speed  of  up  to  a  factor  of  ten 
compared  to  the  CRAY-1  and  a  memory  of  256xl06  64-bit  words.  We  understand  that  it  is  scheduled 
to  be  available  before  the  end  of  1985.  We  also  understand  that  CDC  is  planning  a  successor  to  the 
CYBER  205,  by  1986,  to  consist  of  several  processors,  each  processor  three  to  five  times  as  fast  as  the 
205  processor. 


8 


NRL  REPORT  8820 


There  are  intense  efforts  under  way  to  increase  computing  power  by  two  or  three  orders  of  magni¬ 
tude.  These  research  efforts  are  based  on  the  concept  of  parallel  computing  using  several  different  pro¬ 
cessors.  There  is  an  abundance  of  concepts  for  parallel  computing  systems  as  well  as  special  purpose 
computers  (such  as  finite  element  machines).  We  shall  not  go  into  these  issues  but  instead  refer  to 
Ref.  10  and  the  references  cited  there  for  technical  details.  We  merely  remark  that  this  research  is  of  a 
long-term  nature  and  it  is  not  clear  when  parallel  computing  systems  will  be  available  to  deliver  this 
kind  of  computing  power  for  realistic  problems. 

Let  us  next  try  to  assess  the  feasibility  of  solving  large  problems  efficiently  on  both  the  CRAY-1 
and  the  CRAY-2.  The  following  discussion  is  based  on  estimates  obtained  from  a  simple  scaling  of 
problems  which  have  already  been  solved.  It  is  not  possible  to  predict  how  valid  this  scaling  is. 

For  the  purpose  of  this  discussion,  we  consider  the  Helmholtz  model  described  previously  over  a 
range  of  40  km  at  a  frequency  of  40  Hz.  Furthermore,  we  assume  a  uniform  grid  size  of  eight 
points/wavelength  although,  as  we  have  seen,  we  could  considerably  reduce  the  number  of  equations 
using  nonuniform  grid  sizes.  We  see  from  Table  1  Case  (c)  that  for  a  frequency  of  10  Hz  and  a  range 
of  2  km  the  CPU  time  on  the  VAX  was  2.5  h.  Let  us  now  simply  scale  up  this  problem  to  40  Hz  over 
40  km  multiplying  by  42  for  the  frequency  increase  and  20  for  the  range  increase,  using  the  reasoning 
outlined  above  We  would  then  obtain  for  the  large  problem  (consisting  of  nearly  107  complex  degrees 
of  freedom),  a  CPU  time  of  800  h  on  the  VAX,  20  h  on  the  CRAY-1  and  2  h  on  the  CRAY-2  (assum¬ 
ing  for  the  sake  of  definiteness  that  the  CRAY-2  will  be  a  factor  of  10  faster  than  the  CRAY-1).  This 
rough  scaling  neglects  any  increase  in  iterations  as  the  size  of  the  region  increases. 

There  are  a  number  of  factors  that  cannot  be  predicted  in  advance.  For  one  thing,  it  is  not  really 
possible  to  predict  how  a  given  code  will  behave  on  a  different  computer.  Furthermore,  the  problem  is 
too  large  to  fit  in  the  central  memory  of  the  CRAY-1.  This  necessitates  an  input/output  process  that 
can  cause  a  considerable  increase  in  time.  This  additional  time  can  be  minimized  by  the  use  of  a 
numerical  algorithm  well  designed  to  deal  with  this  difficulty.  This  problem  would  be  considerably  less 
severe  for  the  CRAY-2,  however,  because  of  its  predicted  large  central  memory,  if  we  assume  that  a 
large  part  of  its  memory  is  available  for  this  problem.  Indications  are  that  memory  is  coming  down  in 
price,  so  it  is  to  be  expected  that  a  wide  range  of  computers  will  have  more  memory  available  to  them. 

Finally,  it  is  anticipated  that  the  suggested  improvements  outlined  in  the  previous  section  and, 
more  important,  those  improvements  described  in  the  following  would  lead  to  substantial  improve¬ 
ments  in  the  algorithm  and  computer  code.  It  is  not  possible,  however,  to  predict  the  exact  degree  to 
which  such  improvements  would  reduce  the  storage  requirements  and  CPU  time.  For  these  and  other 
reasons,  the  only  way  adequately  to  assess  the  parameter  ranges  that  can  be  efficiently  modeled  is  by 
means  of  continuing  numerical  studies  employing  the  best  modeling  and  mathematical  techniques  avail¬ 
able. 


Time-Harmonic  Models 

We  now  discuss  mathematical  models  corresponding  to  a  time-harmonic  source.  Such  models 
commonly  occur  in  underwater  acoustic  propagation  problems  [1].  Since  the  harmonic  time  depen¬ 
dence  is  factored  out,  we  are  left  with  an  elliptic  boundary-value  problem.  This  problem  consists  of  a 
system  of  coupled  second-order  elliptic  partial  differential  equations,  each  analogous  to  the  Helmholtz 
equation  considered  previously.  We  proceed  to  discuss  briefly  three  areas  of  mathematical  research  that 
can  have  an  important  bearing  on  the  numerical  solution  of  these  models.  This  is  not  intended  to  be  a 
comprehensive  survey  but  merely  to  indicate  some  of  the  research  directions  that  can  yield  fruitful 
results. 
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Radiation  Boundary  Conditions 


To  solve  the  problem  approximately,  we  must  introduce  an  artificial  vertical  boundary  on  the  right 
side  of  the  computational  domain  (see  Fig.  2)  as  well  as  an  appropriate  radiation  (absorbing)  boundary 
condition  to  simulate  the  outgoing  propagation  of  energy.  This  artificial  boundary  can  intersect  the  fluid 
domain  completely,  the  solid  domain  completely,  or  both  the  fluid  and  the  solid  domains.  Various 
methods  have  been  formulated  for  modeling  radiation  boundary  conditions,  and  much  research  is  con¬ 
tinuing  in  this  area.  The  goal  is  to  formulate  a  boundary  condition  that  minimizes  spurious  reflections 
that  can  contaminate  the  solution  without  increasing  the  size  of  the  computational  domain  unneces¬ 
sarily.  At  the  same  time,  it  is  important  for  this  formulation  to  be  suitable  for  efficient  implementation 
and  solution.  The  choice  of  the  most  efficient  formulation  is  to  a  great  extent  problem  dependent.  For 
our  model  problem  in  the  previous  section,  we  discussed  the  use  of  artificial  dissipation  to  attenuate  the 
wave.  This  method  is  only  applicable  when  no  backscattering  occurs.  Empirical  evidence  indicates  that 
in  such  cases  the  method  can  be  effective  with  only  a  few  wavelengths  of  artificial  dissipation. 

We  first  assume  that  the  artificial  boundary  intersects  only  the  fluid  domain,  so  that  we  consider 
the  acoustic-wave  equation.  In  many  cases,  the  radiation  condition  may  be  expressed  in  terms  of  a 
modal  expansion  for  the  outgoing  solution.  For  example,  for  the  case  of  a  rotationally  symmetric 
medium  the  following  radiation  condition  was  employed  in  combination  with  the  finite  element  method 


H  -  «> 

Here,  the  expression  for  T{u)  involves  a  series  expansion  in  Hankel  functions  to  represent  outgoing 
waves  at  the  artificial  boundary  and  bu/bn  is  the  normal  derivative  at  this  boundary.  This  formulation 
is  valid  whenever  the  problem  is  separable.  This  occurs,  for  example  when  the  sound  speed  is  range 
independent  and  both  the  upper  and  lower  boundaries  are  horizontal.  It  is  only  necessary  to  consider  a 
finite  number  of  terms  in  this  expansion  corresponding  to  the  propagating  modes.  An  analogous 
ary  condition  was  analyzed  in  Ref.  12  in  connection  with  a  variety  of  geometries,  and  it  was  shown  that 
the  finite  element  method  converges  optimally  in  spite  of  the  complicated  nature  of  this  boundary  con¬ 
dition. 

A  hierarchy  of  boundary  conditions  based  on  this  modal  expansion  was  employed  in  Ref.  2.  Each 
successive  boundary  condition  in  this  hierarchy  was  constructed  to  be  exact  for  successively  more  prop¬ 
agation  modes.  Unlike  condition  (4),  these  boundary  conditions  are  local  (i.e.,  the  normal  derivative  of 
w  at  a  point  is  given  in  terms  of  u  and  its  tangential  derivatives  at  that  point).  Hence  the  resulting  sys¬ 
tem  of  equations  is  sparse,  which  leads  to  certain  computational  advantages.  When  there  are  only  one 
or  two  propagating  modes,  these  boundary  conditions  are  more  conveniently  implemented  than  Eq.  (4). 
However,  practical  difficulties  occur  in  the  implementation  of  higher  order  boundary  conditions  due  to 
the  presence  of  higher  order  tangential  derivatives.  The  global  condition  (4)  is  applicable  even  when 
there  are  many  propagating  modes. 

Another  hierarchy  of  local  boundary  conditions  was  developed  and  analyzed  in  Ref.  13  for  the 
time-dependent  wave  equation  using  pseudodifferential  operators.  These  conditions  are  also  applicable 
to  the  Helmholtz  equation  [14],  The  first  boundary  condition  in  this  hierarchy  is  closely  related  to  the 
viscous  boundary  condition,  described  in  Eq.  (5).  It  is  exact  if  the  solution  is  a  plane  wave  normally 
incident  to  the  artificial  boundary.  The  higher  order  boundary  conditions  in  Ref.  13  are  designed  to 
give  better  approximations  as  the  deviation  from  normal  incidence  increases.  As  before,  however, 
there  are  computational  difficulties  in  implementing  the  high-order  boundary  conditions,  since  they 
involve  high-order  derivatives. 

Next,  suppose  that  the  artificial  boundary  intersects  the  solid  domain.  A  convenient  method  for 
handling  the  absorbing  boundary  condition  in  conjunction  with  the  displacement  formulation  of  the 
elliptic  wave  equation  was  described  in  Ref.  15.  A  plane  wave  boundary  condition  of  the  form 
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<r„  =  pcdv„,  <r,  =  pcsv,  VJ' 

is  applied  at  the  artificial  boundary.  Here,  <r„  and  cr,  are  the  normal  and  tangential  interface  stresses,  p 
is  the  solid-medium  mass  density,  cd  and  cs  are  the  dilatational  and  shear  wave  speeds,  and  v„  an<^  v» 
are  the  corresponding  normal  and  tangential  velocities  on  the  solid  boundary  surface.  Conditions  ( .  ) 
are  exact  if  the  radiating  energy  at  the  surface  consists  of  plane  waves  that  are  normally  incident  to  the 
surface.  It  was  pointed  out  in  Ref.  15  that  these  boundary  conditions  can  still  be  accurate  even  when 
this  condition  does  not  hold.  The  artificial  boundary  acts  as  if  viscous  dampers  have  been  applied  nor¬ 
mally  and  tangentially  to  the  solid-domain  boundary  points.  When  the  dominant  waves  traveling 
through  the  media  were  Rayleigh  waves  rather  than  plane  waves,  an  analogous  boundary  condition  was 
applied  in  Refs.  15  and  16.  The  boundary  conditions  in  Eq.  (5)  were  also  applied  to  treat  the  case  in 
which  the  artificial  boundary  intersects  both  the  fluid  and  solid  media  [9].  Empirical  evidence  indicates 
that  these  boundary  conditions  can  give  satisfactory  results  if  they  are  located  far  enough  from  the 
source  force  functions  [17,18]. 


Finally,  we  note  that  there  are  various  alternative  methods  for  treating  the  unbounded  region, 
such  as  the  use  of  infinite  elements  [19].  In  this  case,  the  usual  piecewise  polynomial  basis  functions  are 
replaced  by  special  functions  in  the  outermost  layer  of  elements.  These  infinite  elements  are  chosen  so 
as  to  simulate  the  behavior  of  the  solution  near  infinity.  Another  approach  for  treating  unbounded 
domains  consists  of  coupling  the  finite  element  method  with  an  integral  operator  on  the  outer  boundary 
[20].  The  Green’s  function  for  the  integral  operator  simulates  the  behavior  of  the  solution  near 
infinity.  For  additional  methods  for  treating  radiation  boundary  conditions,  see  Refs.  9  and  14. 


Adaptive  Discretization  Methods 

To  solve  the  resulting  elliptic  boundary-value  problem  approximately,  we  must  discretize  this 
problem.  Generally,  this  may  be  efficiently  done  by  a  variety  of  finite  difference  or  finite  element 
methods.  Finite  element  methods  are  better  suited  to  treating  the  complicated  boundaries  and  bound¬ 
ary  conditions  that  often  arise  in  the  propagation  problems  considered  here.  Because  of  their  great 
flexibility  in  modeling  complicated  problems,  we  recommend  the  use  of  finite  element  methods.  It  is 
not  clear,  at  the  present  time,  whether  high-order  or  low-order  methods  would  be  more  efficient  for 
these  problems.  This  question  needs  to  be  investigated  further. 

Adaptive  computational  methods  have  recently  proved  to  be  a  powerful  tool  in  the  numerical 
solution  of  partial  differential  equations,  and  considerable  research  is  continuing  along  these  lines. 
These  algorithms  have  built  into  them  convenient  methods  for  obtaining  measures  of  accuracy  and 
adapting  the  discretization  automatically  to  the  evolving  solution,  so  as  to  obtain  a  desired  level  of  accu¬ 
racy  with  a  minimum  of  arithmetic  operations.  These  algorithms  automatically  provide  for  smaller 
mesh  sizes  in  a  region  where  the  solution  is  not  smooth  (e.g.,  where  there  is  a  large  velocity  or  density 
gradient)  and  larger  mesh  sizes  when  the  solution  is  smooth.  Recently  developed  adaptive  discretiza¬ 
tion  methods  have  resulted  in  significant  improvements  in  the  numerical  solution  of  elliptic  boundary- 
value  problems  in  elasticity  and  other  areas  (see  Refs.  21  and  22,  and  the  references  cited  there.) 

We  feel  that  the  use  of  adaptive  methods  can  also  lead  to  dramatic  improvement  in  the  numerical 
solution  of  underwater  acoustic  propagation  problems.  As  we  have  demonstrated  (in  a  previous  sec¬ 
tion,  in  Ref.  6,  and  in  Ref.  23),  the  use  of  nonuniform  grid  sizes  in  connection  with  the  Helmholtz 
equation  can  lead  to  improved  accuracy  with  considerably  fewer  grid  points.  Adaptive  methods  provide 
an  efficient  means  of  changing  the  mesh  sizes  (and  other  important  parameters)  systematically  in  a 
nearly  optimal  manner.  Furthermore,  the  adaptive  method  provides  a  means  for  assessing  the  accuracy 
of  the  computed  solution.  Hence,  another  potential  advantage  is  that  the  error  due  to  various 
simplified  propagation  models  may  be  conveniently  determined. 
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Iterative  Methods 

The  finite  element  or  finite  difference  discretization  of  our  model  results  in  a  sparse,  indefinite, 
non-self-adjoint  system  of  linear  equations.  As  we  have  mentioned  previously,  iterative  solution 
methods  are  far  superior  to  direct  methods  for  very  large  problems.  However,  standard  iterative 
methods  are  not  applicable  to  systems  that  are  indefinite  and  non-self-adjoint.  The  iterative  method 
discussed  previously  in  connection  with  the  Helmholtz  equation  is  also  applicable  to  the  more-general 
boundary-value  problems  considered  here.  It  is  based  on  the  conjugate  gradient  method  applied  to  the 
normal  equations  with  an  appropriate  preconditioning  (see  Ref.  2).  The  performance  of  the  iterative 
method  has  a  significant  bearing  on  the  cost  of  the  computation.  Hence,  the  study  of  preconditioned 
conjugate  gradient  methods  and  other  iterative  methods  in  connection  with  elliptic  boundary-value 
problems  is  receiving  a  great  deal  of  attention  (see,  e.g„  Refs.  24  and  25).  We  shall  briefly  consider 
ways  of  potentially  improving  the  iterative  method  described  previously. 

The  choice  of  preconditioner  in  the  preconditioned  coryugate  gradient  method  is  crucial  in  deter¬ 
mining  the  number  of  iterations  required  for  a  desired  accuracy  and,  hence,  the  amount  of  computa¬ 
tion.  Another  important  factor  is  the  method  in  which  the  preconditioner  is  implemented  in  the  conju¬ 
gate  gradient  method.  Thus  far  we  have  only  implemented  a  preconditioner  based  on  one  sweep  of 
SSOR  in  connection  with  the  problem  of  Eqs.  (1)  and  (2).  As  demonstrated  in  Refs.  2  and  26,  the 
number  of  iterations  using  this  preconditioner  grows  as  O  (A-1)  as  the  mesh  size  h  — 1 ►  0  There  are 
several  alternative  preconditioners  that  have  proved  to  be  quite  effective  in  connection  with  the  conju¬ 
gate  gradient  method  [24],  Some  of  these  are  currently  being  investigated  in  connection  with  the 
Helmholtz  equation  [26],  and  in  theory  appear  to  exhibit  a  faster  convergence  rate  than  SSOR  for  small 
mesh  sizes.  For  example,  it  can  be  seen  theoretically  that  a  preconditioner  based  on  a  fast  solver  (such 
as  a  multigrid  method)  has  a  convergence  rate  that  is  independent  of  h.  Hence  the  number  of 

iterations  is  indpendent  of  h  and  such  a  preconditioner  would  appear  to  be  superior  for  sufficiently  small 
mesh  sizes. 

It  remains  to  be  seen  whether  preconditioners  that  are  more  effective  as  h  —  0  will  be  more 
effective  in  realistic  underwater  acoustic  problems.  There  are  many  other  factors  that  play  a  role  in  the 
choice  of  preconditioners,  such  as  their  behavior  with  respect  to  the  wave  number  and  boundary  condi¬ 
tions.  Finally,  we  recall  that,  to  utilize  fully  the  capabilities  of  vector  computers  such  as  the  CRAY-1 
or  CYBER  205,  the  algorithm  should  be  vectorizable.  This  is  not  the  case  of  the  preconditioners  we 
have  implemented  thus  far.  For  examples  of  preconditioned  conjugate  gradient  methods  that  have 
been  vectorized,  see  Refs.  27  and  28.  The  choice  of  an  optimal  preconditioner  for  a  given  class  of 

problems  and  a  given  computer  can  be  complicated,  but  it  can  be  important  in  reducing  the  cost  of  the 
computation. 

Time-Dependent  Models 

Time-dependent  models  are  most  appropriate  when  the  structure  is  subjected  to  a  transient  force 
as  a  pulse.  This  is  commonly  the  case,  for  example,  in  seismology.  The  mathematical  model  in  such 
cases  is  an  initial-boundary-value  problem  for  a  coupled  system  of  wave  equations.  This  hyperbolic 
problem  is  different  mathematically  from  the  elliptic  problem  associated  with  the  time-harmonic  model 
and  hence  requires  different  computational  methods  for  its  solution.  We  shall  very  briefly  discuss  some 
of  these  methods  as  well  as  some  ways  of  potentially  speeding  up  the  computations. 

Boundary  Conditions 

As  for  the  time-harmonic  model,  there  are  different  formulations  of  the  radiation  condition  on  an 
artificial  outer  boundary.  It  is  important  to  choose  the  most  appropriate  one  for  a  given  problem.  Typi¬ 
cal  formulations  for  the  time-dependent  model  are  based  on  local,  approximate  absorbing  boundary 
conditions.  Viscous  damper  boundary  conditions  such  as  those  described  in  Eq.  (5),  are  commonly 
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used.  However,  they  require  that  the  artificial  boundary  be  placed  sufficiently  far  away  from  the  source 
of  excitation  so  that  reflections  are  not  pernicious.  A  hierarchy  of  local,  approximate  boundary  condi¬ 
tions  was  developed  in  Ref.  13  for  the  acoustic-wave  equation  and  in  Ref.  29  for  the  elastic- wave  equa¬ 
tion.  As  in  the  case  of  the  viscous  boundary  condition,  the  accuracy  of  these  boundary  conditions 
depends  on  the  deviation  of  the  wave  from  normal  incidence  at  the  artificial  boundary.  The  higher 
order  conditions  are  more  accurate,  and  hence  the  size  of  the  computational  domain  may  be  decreased. 
However,  they  are  more  difficult  to  implement. 

Finally,  we  mention  a  boundary-condition  approach  developed  in  Ref.  30  to  eliminate  reflections 
for  transient  problems.  This  method  is  based  on  the  calculation  of  two  independent  solutions  in  which 
the  reflections  are  of  opposite  sign.  The  addition  of  these  solutions  cancels  the  reflections,  leaving  only 
the  energy  originally  incident  on  the  boundary.  The  two  solutions  may  be  obtained,  for  example,  by 
appling  homogeneous  Dirichlet  and  Neumann  boundary  conditions. 

Discretization  Methods 

We  next  consider  the  question  of  discretizing  the  resulting  initial-boundary-value  problem.  Now 
we  must  consider  the  time  discretization  in  addition  to  the  spatial  discretization.  As  before  the  spatial 
variables  may  be  discretized  using  either  a  finite  difference  or  finite  element  method.  As  for  the  time 
discretization,  empirical  experience  indicates  that  explicit  finite  difference  time  discretizations  are  more 
efficient  than  implicit  ones  in  spite  of  stability  constraints  (see  Refs.  8  and  31),  since  the  latter  method 
involves  the  solution  of  linear  equations  at  each  time  step.  Explicit  formulations  express  the  value  of  a 
variable  at  some  point  at  a  future  time  in  terms  of  the  value  of  the  variable  at  that  point  and  neighbor¬ 
ing  points  at  the  present  time  and  past  times.  Hence,  it  is  no  longer  necessary  to  solve  large  systems  of 
linear  equations.  The  large  computation  time  for  these  problems  is  due  mainly  to  the  large  number  of 
time  steps  required  and  the  large  number  of  spatial  variables. 

We  propose  two  methods  for  reducing  these  computational  costs.  The  first  suggestion  is  the  use 
of  an  adaptive  discretization  method.  In  a  previous  section  we  described  the  potential  advantages  asso¬ 
ciated  with  such  methods.  Adaptive  methods  appropriate  for  hyperbolic  problems  are  different  from 
those  most  useful  for  elliptic  problems.  See  Ref.  32  and  the  references  cited  there  for  a  discussion  of 
adaptive  discretization  methods  for  time-dependent  problems.  The  second  suggestion  is  the  use  of 
higher  order  discretization  methods  in  the  spatial  variables.  As  demonstrated  in  Ref.  33,  significant 
gains  can  be  expected  from  the  use  of  higher  order  methods.  These  include  higher  order  finite 
difference  and  finite  element  methods  as  well  as  spectral  methods.  Recent  work  [34]  dealing  with  the 
discretization  error  as  the  number  of  wavelengths  increases  (i.e.,  at  high  frequencies  or  in  large 
domains)  indicates  that,  particularly  in  these  cases,  higher  order  spatial  discretization  can  yield  substan¬ 
tial  improvements  for  both  transient  and  stationary  propagation  models.  On  the  other  hand,  as  indi¬ 
cated  in  Ref.  33,  it  often  suffices  to  use  a  second-order  time-discretization  method. 

CONCLUSIONS 

We  have  considered  the  numerical  solution  of  direct  deterministic  underwater  acoustic- 
propagation  problems  using  finite  difference  and  finite  element  methods.  We  first  discussed  results 
obtained  after  implementing  a  finite  element  code  on  a  VAX  11/780  and  modifying  it  to  treat  a  propa¬ 
gation  model  based  on  the  two-dimensional  Helmholtz  equation  with  a  variable  sound  speed.  An 
important  feature  of  the  code  is  the  implementation  of  a  recently  developed  iterative  method  based  on 
the  preconditioned  conjugate  gradient  method  for  solving  the  resulting  large,  sparse,  indefinite,  non- 
self-adjoint  system  of  linear  equations.  Since  large  matrices  do  not  have  to  be  stored  or  inverted,  we 
were  able  to  solve  efficiently  over  35,000  complex  equations.  Furthermore,  simple  changes  in  the  code 
were  outlined  that  would  substantially  increase  both  storage  capabilities  and  speed. 
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We  next  discussed  various  issues  involved  in  the  numerical  solution  of  general  two-  and  three- 
dimensional  propagation  models.  Such  models  involve  the  coupling  of  the  acoustic-wave  equation  in  a 
fluid  with  the  elastic-wave  equation  in  a  solid  and  the  specification  of  appropriate  interface  and  bound¬ 
ary  conditions.  This  allows  for  the  computational  modeling  of  general  ocean  environments,  including 
complicated  bottom  structures  as  well  as  a  layer  of  ice  on  the  surface.  We  considered  time-harmonic 
(stationary)  formulations,  which  commonly  occur  in  underwater  acoustic  models,  as  well  as  time- 
dependent  formulations  corresponding  to  a  transient  force  such  as  a  pulse.  The  resulting  models  may 
be  solved  approximately  by  finite  difference  or  finite  element  methods. 

The  numerical  methods  for  treating  stationary  and  transient  problems  are  very  different.  For  tran¬ 
sient  problems,  explicit  finite  difference  time  discretizations  appear  to  be  more  efficient  than  implicit 
ones,  since  the  latter  involve  the  solution  of  linear  equations  at  each  time  step.  For  stationary  prob¬ 
lems,  we  recommend  the  use  of  iterative  methods  over  direct  methods.  For  both  stationary  and  tran¬ 
sient  formulations  we  recommend  the  use  of  finite  element  methods  for  the  spatial  discretization 
because  of  their  great  flexibility  in  modeling  complicated  problems. 

The  size  of  the  problem  depends  to  a  large  extent  on  the  frequency  and  the  size  of  the  computa¬ 
tional  domain.  In  particular,  because  the  waves  must  be  resolved  in  all  directions,  the  number  of 
degrees  of  freedom  in  two  (or  three)  dimensions  generally  increases  as  at  least  the  square  (or  cube)  of 
the  frequency.  The  size  of  the  problem  that  can  be  efficiently  treated  depends  on  the  computer  power 
available  as  well  as  the  mathematical  modeling  techniques  and  the  numerical  algorithm  employed.  Vec¬ 
tor  computers  are  expected  to  improve  current  computer  power  by  approximately  a  factor  of  ten  over 
the  next  few  years.  Longer  term  research  in  multiprocessor  systems  is  intended  to  result  in  improve¬ 
ments  of  two  and  three  orders  of  magnitude.  The  numerical  methods  we  have  discussed  based  on 
finite  elements  or  finite  differences  are  well  suited  to  parallel  computation. 

Furthermore,  we  described  several  mathematical  modeling  and  numerical  investigations  that  are 
expected  to  result  in  substantial  improvements  in  computational  efficiency.  These  include  the  use  of 
methods  for  modeling  the  radiation  boundary  condition  so  as  to  reduce  the  size  of  the  truncated 
computational  domain,  the  use  of  techniques  for  improving  the  iterative  method,  and  the  use  of  adap¬ 
tive  discretization  methods.  Numerical  methods  can  also  be  useful  for  assessing  the  accuracy  of  various 
simplified  propagation  models.  We  feel  that,  at  the  present  time,  a  variety  of  complicated,  realistic, 
two-dimensional  propagation  problems  can  be  efficiently  solved  by  numerical  methods  without  the 
necessity  of  resorting  to  simplified  models.  Furthermore,  in  view  of  anticipated  improvements  in  com¬ 
puter  power  and  numerical  methods,  it  is  likely  that  the  same  will  be  true  for  three-dimensional  and 
much  larger  two-dimensional  problems  within  a  few  years. 
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